#=
观测量
=#


#=
d lnZ / dp = d Z / dp / Z
= d ∫dR ρ(R) / dp / ∫dR ρ(R)
= ∫dR ρ(R) (1/ρ(R)) dρ(R)/dp / ∫dR ρ(R)
= ∫dR ρ(R) dlnρ(R)/dp / ∫dR ρ(R)
= { ∫dR |ρ(R)| ρ(R)/|ρ(R)| O(R) / ∫dR |ρ(R)| } / { ∫dR |ρ(R)| ρ(R)/|ρ(R)| / ∫dR |ρ(R)| }
=#

#=
Z = ∑_{P} sgn(P) ∫dr1 <r1| e^{-βH} |Pr1>
= ∫dr1dr2...drm <r1| e^{-ΔβH} |r2> <r2| e^{-ΔβH} |r3> ... <rm| e^{-ΔβH} |Pr1>
= ∫dR ρ(r1,r2) ρ(r2,r3) ... ∑_{P} sgn(P) <rm| e^{-ΔβH} |Pr1>
= ∫dR ρ(r1,r2) ρ(r2,r3) ... ∑_{P} sgn(P) e^{-ΔβU(rm,Pr1)} e^{-π(rm-Pr1)^2/λe^2}
##势能要对所有的pair求和，所以不依赖交换 U(rm,Pr1) = U(rm,r1)
= ∫dR ρ(r1,r2) ρ(r2,r3) ... e^{-ΔβU(rm,r1)} ∑_{P} sgn(P) e^{-π(rm-Pr1)^2/λe^2}
= ∫dR ρ(r1,r2) ρ(r2,r3) ... e^{-ΔβU(rm,r1)} det|Ψ|
= ∫dR ρ(r1,r2) ρ(r2,r3) ... ρ̃(r_m,r_1)
=#

#=
dlnρ(R) = dln { ∏^{M-1}_1 ρ(r_i,r_{i+1}) } ρ̃(r_m,r_1)
= ∑^{M-1}_1 dln ρ(r_i,r_{i+1}) + dln ρ̃(r_m,r_1)
= ∑^{M-1}_1 dln e^{-ΔβU(r_i,r_{i+1})}e^{-π(r_{i+1}-r_i)^2/λe^2} + dln e^{-ΔβU(r_m,r_1)}det|Ψ|
= ∑^{M-1}_1 d -ΔβU(r_i,r_{i+1}) -π(r_{i+1}-r_i)^2/λe^2 + d -ΔβU(r_m,r_1) + dln det|Ψ|
=#

#=
1) dln det|Ψ| = tr |Ψ|^-1 d|Ψ|
2) dln det|Ψ| / dp = ∑_{i,j} (d Ψij / dp) { dln det|Ψ| / dΨij }
= ∑_{i,j} (d Ψij / dp) { tr |Ψ|^-1 d|Ψ| / dΨij }
= ∑_{i,j} (d Ψij / dp) { ∑_{x,y} {|Ψ|^-1}xy d Ψyx / dΨij }
= ∑_{i,j} (d Ψij / dp) {|Ψ|^-1}ji
= tr {d|Ψ|/dp} {|Ψ|^-1}
##这里的逆矩阵={余子式的转置}/det|Ψ|
## 另一种求导方式：再每一列上求导
## Ψij dΨij/dp Ψij
## Ψij dΨij/dp Ψij
## Ψij dΨij/dp Ψij
## det|Ψ|^-1 d ∑sgn(P) ψ...dΨ/dp...
## 和上面的逆矩阵*dΨ/dp结果是一样的
=#

#=
βE = -βdlnρ(R)/dβ
= -β{ d( ∑^m_1 d -ΔβU(r_i,r_{i+1}) + ∑^{m-1}_1 -π(r_{i+1}-r_i)^2/λe^2) / dβ + {d|Ψ|/dβ} {|Ψ|^-1}^T }
= -MΔβ d ∑-ΔβU /dMΔβ - MΔβ d-πδr^2/λe^2 /dMΔβ  - MΔβ tr {d|Ψ|/dMΔβ} {|Ψ|^-1}^T 
= Δβ d ∑ΔβU /dΔβ + Δβ dπδr^2/λe^2 /dΔβ  - Δβ tr {d|Ψ|/dΔβ} {|Ψ|^-1}^T 
=#

#=
dΨij/dβ = de^{-π(ri-rj)^2/λe^2}/dβ
## 如果利用ri_m = ri_{1} + yi λe
## = de^{-π(ri-rj+yi λe)^2/λe^2}/dβ
## = de^{-π∑_{x,y,z} (ri-rj)_x^2/λe^2 + (yi^2λe^2)_x/λe^2 + 2(ri-rj*yλe)_x/λe^2}/dβ
## 第二个没有了
## = d-π∑_{x,y,z} {(ri-rj)_x^2/λe^2 + 2(ri-rj*yλe)_x/λe^2}/MdΔβ * e^{-π(ri-rj+yi λe)^2/λe^2}
## = d-π∑_{x,y,z} {(ri-rj)_x^2 m/2πh̄^2Δβ + 2(ri-rj*y)_x sqrt(m/2πh̄^2Δβ)}/MdΔβ * e^{-π(ri-rj+yi λe)^2/λe^2}
## = -π∑_{x,y,z} {-(ri-rj)_x^2 m/2πh̄^2Δβ^2 - (|ri-rj|*y)_x sqrt(m/2πh̄^2Δβ)/Δβ}/M * e^{-π(ri-rj+yi λe)^2/λe^2}
## = -π∑_{x,y,z} {-(ri-rj)_x^2 /λe^2 - 2(|ri-rj|*yλe)_x /λe^2}/MΔβ * e^{-π(ri-rj+yi λe)^2/λe^2}
## = +π(ri-rj)⋅(ri-rj+y)/λe^2/MΔβ * e^{-π(ri-rj+yi λe)^2/λe^2}
# 上面y和yλe没做区分，实际上应该用y=ξλe，最后还是都吸收回来变成y了
``````
如果ri-rj+yλe触发了周期边界条件变成ri-rj+yλe±L,则让ri-rj变成ri-rj±L,其余的不变
=#

#=
dΨij/dβ = de^{-π(ri-ri+1)^2/λe^2}/dβ
## 如果ri和ri+1之间的距离用λe为单位，d -πy^2/ dβ = 0
=#

"""
对ψ的矩阵元进行求导(对Δβ求导并乘以Δβ)
rpos就是最后一个时刻的位置，r0是第一个
"""
function exch_derv_beta(rpos::Matrix{T}, r0::Matrix{T}, vol::Vector{Float64};
    λe::Float64=_ThermWL) where {T}
    Np, Nd = size(rpos)
    dervmat = zeros(T, Np, Np)
    for i in Base.OneTo(Np); for j in Base.OneTo(Np)
        if i == j
            #对角的元素没有导数，(ri-rj)⋅(ri-rj+y) = 0
            continue
        end
        rimrj = r0[i, :] - r0[j, :]
        rimrjpy = rpos[i, :] - r0[j, :]
        #
        disij = zeros(Nd)
        for id in Base.OneTo(Nd)
            disij[id] = dis_to_1st(rimrjpy[id], vol[id])
        end
        rimrjpy = rimrjpy - disij
        rimrj = rimrj - disij
        #
        rinn = sum(rimrj .* rimrjpy)
        rsqu = sum(rimrjpy.^2)
        dervmat[i, j] = π*rinn/λe/λe*exp(-π*rsqu/λe/λe)
    end; end
    return dervmat
end


#=
Z = ∫dR ρ(R)

dlnZ / dp = ∫dR lnρ(R) / dp

lnρ(R) = dln { ∏^{M-1}_1 ρ(r_i,r_{i+1}) } ρ̃(r_m,r_1)
=∑^{M-1}_1 d -ΔβU(r_i,r_{i+1}) -π(r_{i+1}-r_i)^2/λe^2 + d -ΔβU(r_m,r_1) + dln det|Ψ|
##后面的交换项在前面处理，第二个动能项没有贡献
#λ11 = sqrt((_hbar^2)*_Δβ/_ElecMass)
#Kelbg
#( 1-e^{-|r|^2/λ12^2} + √π |r|/λ12 (1-erf(|r|/λ12)) )/|r|
#对Kelbg的导数是
(1/|ri-rj| - <ri-rj|yi-yj>/2|ri-rj|^3)(1-e^{-|ri-rj|^2/λ11^2}) + √π/2λ11 (1-erf(|ri-rj|/λ11))
=#

"""
Kelbg随beta的导数(对Δβ求导并乘以Δβ)
被求导的对象应该是ΔβΦee，因为修改了Δβ以后β也会随着修改
"""
function _Kelbg_derv_beta(rij::Vector{T}, r0ij::Vector{T}, λab::Float64) where {T}
    absrij = sqrt(sum(rij.^2))
    #ri = r0i + yi (yi=ξiλe)
    yij = rij - r0ij
    rijcyij = sum(rij .* yij)
    #(1/|ri-rj| - <ri-rj|yi-yj>/2|ri-rj|^3)
    #c1 = rijcyij/2/(absrij^3)#1/absrij - rijcyij/2/(absrij^3)
    c1 = 1/absrij - rijcyij/2/(absrij^3)
    #(1-e^{-|ri-rj|^2/λ11^2})
    c2 = 1 - exp(-(absrij^2)/(λab^2))
    #√π/2λ11 (1-erf(|ri-rj|/λ11))
    c3 = (√π)/2/λab * (1 - erf(absrij/λab))
    #c3 = (√π)/λ11 * (1 - erf(absrij/λ11))
    #
    if isapprox(0.0, c2, atol=eps())
        return c3
    else
        return c1*c2 + c3
    end
    #λ12 = sqrt((_hbar^2)*_Δβ/_ElecMass)
    #(1 -e^{-|r|^2/λ12^2} +√π |r|/λ12 (1-erf(|r|/λ12)) )/|r|
    #第一项对Δβ求导是0
    #第二项 -e^{-_ElecMass|r|^2/(_hbar^2)*_Δβ}
    # -e^{-_ElecMass|r|^2/(_hbar^2)*_Δβ} ∂(-_ElecMass|r|^2/(_hbar^2)*_Δβ)/∂Δβ
    # -e^{-_ElecMass|r|^2/(_hbar^2)*_Δβ} _ElecMass|r|^2/(_hbar^2)*Δβ/Δβ
    # -e^{-|r|^2/λ12^2} |r|^2 / λ12^2 / Δβ
    #第三项√π |r|/λ12
    #d √π |r| sqrt(_ElecMass) / sqrt((_hbar^2)*_Δβ) / dΔβ
    # - ( √π |r| sqrt(_ElecMass) / (_hbar^2)*_Δβ ) ( d sqrt((_hbar^2)*_Δβ) / dΔβ )
    # - ( √π |r| sqrt(_ElecMass) / (_hbar^2)*_Δβ ) _hbar / 2sqrt(Δβ)
    # - ( √π |r| sqrt(_ElecMass) / _hbar*sqrt(Δβ) ) / 2*Δβ
    # - ( √π |r|/λ12 ) / 2*Δβ
    #第四项
    #-√π |r|/λ12 erf(|r|/λ12))
    #( √π |r|/λ12 ) erf(|r|/λ12)) / 2*Δβ + -√π |r|/λ12 2/√π e^{-|r|^2/λ12^2} (d |r|/λ12 / dΔβ)
    #( √π |r|/λ12 ) erf(|r|/λ12)) / 2*Δβ - 2|r|/λ12 e^{-|r|^2/λ12^2} *( -|r|/λ12 /2*Δβ )
    #( √π |r|/λ12 ) erf(|r|/λ12)) / 2*Δβ + |r|^2/λ12^2 e^{-|r|/λ12} / Δβ
    #
    #都合并起来以后
    #-e^{-|r|^2/λ12^2} |r|^2 / λ12^2 / Δβ就抵消了
end


"""排斥"""
function Kelbg_repulse_dbeta(rij::Vector{T}, r0ij::Vector{T}; λ11=_WLElec) where {T}
    return _Kelbg_derv_beta(rij, r0ij, λ11)
end


"""吸引"""
function Kelbg_attract_dbeta(rij::Vector{T}, r0ij::Vector{T}; λ22=_WLNeur, Cne=_NeurCoul) where {T}
    return -Cne*_Kelbg_derv_beta(rij, r0ij, λ22)
end

